%问题二代码

global W_a ;
global W_b ;
global h2;
global h1;
global L0;
%W_a = 6;
%W_b = 6;

%在此处是什么意思？
WAAA=5.5:0.1:6.5;
WBBB=5.5:0.1:6.5;%起始值是 5.5，步长是 0.1，结束值是 6.5

R_year=[];
E_year=[];

for W_a = WAAA
    for W_b = WBBB
        if W_a<W_b
            continue;
        else
            h2_range = linspace(W_b/2,6,5);%生成 2.75 到 6 生成 5 个点
            %h2_range =2.7500    3.5500    4.3500    5.1500    5.9500
            for h2 = h2_range
                h1 = 80;
                L0_range=linspace(0,350,5);
                for L0 = L0_range
                    %入射角
                    %%计算太阳赤纬角
                    D=[-59 -28 0 31 61 92 122 153 184 214 245 275];
                    R_total = [];
                    EEEE_total=[];
                    EEEE=[];
                    for i = 1:1:size(D,1)
                        %disp(i);
                        [R,EEEE]=Mai(D(i));%将日期传进去
                        %得到了新的向量
                        R_total=[R_total;R];
                        EEEE_total=[EEEE_total;EEEE];
                    end
                    R_year=[R_year;mean(R_total)];%计算平均值
                    E_year=[E_year;mean(EEEE_total) h2 L0];
                end
            end
        end
    end
end
date = data_func;
function[R,EEEE] =  Mai(Day)
global W_a;
global W_b;
global h2;
global h1;
global L0;
EEEE=[];
D=Day;
STT=[9 10.5 12 13.5 15];
result=[];%日结果详细矩阵5*4个时间点
eta_ref=0.92;%镜面反射率
O_R=[0 0 h1];%中心坐标；相当于是永远是以它自身作为坐标原点的
data = readtable("附件.xlsx");%定日镜坐标
data=table2array(data);
%计算太阳赤纬角
sin_delta=sin(2*pi*D/365)*sin(2*pi/360*23.45);
delta=asin(sin_delta);
phi=deg2rad(39.4);
for XH_ST=1:size(STT,2)
    ST=STT(XH_ST);
    omega=pi/12*(ST-12);%太阳时角
    %计算sin(alpha_s)
    sin_alpha_s = cos(delta) * cos(phi) * cos(omega) + sin(delta) * sin(phi);
G_0=1.366; %太阳常数
Hh=3;%3km
ha=0.4237 - 0.00821*(6 - Hh)^2;%大气压
hb=0.5055 + 0.00595*(6.5 - Hh)^2;%大气压
hc=0.2711 + 0.01858*(2.5 - Hh)^2;%大气压
DNI = G_0*(ha + hb*exp(-hc/sin_alpha_s)); %直接法辐照度
alpha_s = asin(sin_alpha_s); % 得到 alpha_s 的弧度值
% 根据公式计算 cos(gamma_s)
cos_gamma_s = (sin(delta) - sin_alpha_s * sin(phi)) / (cos(alpha_s) * cos(phi));
gamma_s = acos(cos_gamma_s); % 得到 gamma_s 的弧度值
x=cos(alpha_s)*sin(gamma_s);
y=cos(alpha_s)*cos(gamma_s);
z=sin(alpha_s);
S_i=[x y z];%太阳入射方向单位向量的相反
Store=zeros(size(data,1),4);%光学效率 阴影遮挡效率 余弦效率 截断效率
E_field=0;%定日镜场的输出热功率
for num=1:1:size(data,1)%%定日镜编号
disp(num);
O_A=[data(num,1) data(num,2) h2];%定日镜中心坐标
%计算集热器截断效率 eta_trunc
d_HR=sqrt((O_A(1)-O_R(1))^2+(O_A(2)-O_R(2))^2+(O_A(3)-O_R(3))^2);%定日镜镜面中心与集热器中心的距离
eta_at = 0.99321 - 0.0001176*d_HR + 1.97e-8 * (d_HR^2); % 计算吸收塔的吸收率
S_AR=(O_R-O_A)/norm(O_R-O_A);%定日镜中心到集热器中心的单位向量 %% 反射角
% 计算法向量
S_n=(S_i+S_AR)/norm(S_i+S_AR);%法向量单位向量
% 计算方位角 Azimuth (θ)
A_H = atan2(S_n(1), S_n(2));
% 计算俯仰角 Elevation (φ)
E_H = acos(S_n(3));
M1 = [1, 0, 0; 0, cos(E_H), sin(E_H); 0, -sin(E_H), cos(E_H)];
M2 = [cos(pi-A_H), sin(pi-A_H), 0; -sin(pi-A_H), cos(pi-A_H), 0; 0, 0, 1];
% 计算矩阵乘积
T = M1 * M2; %镜坐标系变为地坐标系的转换矩阵
S_An=[0 0 1];%镜面 A 坐标系上的的法向量
SSS=S_An*T;
tolerance = 1e-10;
difference = abs(SSS - S_n);
if all(difference(:) >tolerance)
error('法向量转换错误');
end
%
%计算余弦效率 eta_cos
eta_cos = dot(S_i, S_n) / (norm(S_i) * norm(S_n)); % 计算入射向量相反和法向量之间的夹角的余弦值
[P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =shadow_range(-S_i, h2) ;%i 是-S_i
isInside = pointInRectangle(O_A(1), O_A(2), P_top_1_intersect, P_top_2_intersect,P_bot_1_intersect, P_bot_2_intersect);
if isInside %在里面
eta_sb=0;
eta_trunc=0;
Store(num,1)=0;
Store(num,2)=0;
Store(num,3)=eta_cos;
Store(num,4)=0;
else %不在里面
R = 100;
X = linspace(-1/2 * W_a,1/2 * W_a, 5);
Y = linspace(-1/2 * W_b,1/2 *W_b, 5);
% Assuming S_i, S_AR, and T are defined elsewhere
T_inv = inv(T);
% 在指定的半径内寻找定日镜
distances = sqrt((data(1:num-1, 1) - O_A(1)).^2 + (data(1:num-1, 2) - O_A(2)).^2);
valid_indices = find(distances <= R);
% Main loops
valid_point_count = 0;%有效点
valid_points = 0;%每个点的有效光线
obstructed_mirrors = [];
for x = X
for y = Y
dot_A = [x, y, 0];
dot_AG = dot_A * T + O_A;
is_valid = true;
for idx = valid_indices'
O_B = [data(idx, 1), data(idx, 2), h2];
dot_B = (dot_AG - O_B)*T_inv;
S_iB = S_i*T_inv;
H2 = [(S_iB(3)*dot_B(1) - S_iB(1)*dot_B(3))/S_iB(3),(S_iB(3)*dot_B(2) - S_iB(2)*dot_B(3))/S_iB(3), 0];
if ~(-1/2*W_a <= H2(1) && H2(1) <= 1/2*W_a && -1/2*W_b <= H2(2) && H2(2)<= 1/2*W_b)%不在该平面内为真
S_ARB = S_AR*T_inv;
H1 = [(S_ARB(3)*dot_B(1) - S_ARB(1)*dot_B(3))/S_ARB(3),(S_ARB(3)*dot_B(2) - S_ARB(2)*dot_B(3))/S_ARB(3), 0];
if -1/2*W_a <= H1(1) && H1(1) <= 1/2*W_a && -1/2*W_b <= H1(2) &&H1(2) <= 1/2*W_b %在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B,distances(idx)];
is_valid = false;
break;
end
else%在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B, distances(idx)];is_valid = false;
break;
end
end
if is_valid
valid_point_count = valid_point_count + 1;
%这个点是有效的计算集热器接收到的光
% 计算入射向量和法向量之间的夹角
cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
S_nt = S_n / cosTheta_i;
d_ki = S_nt - S_AR;
d_ki = d_ki / norm(d_ki);
theta_range = linspace(-4.65e-3,4.65e-3, 5);
theta2_range = linspace(0, 2*pi, 12);
% 循环 theta 值
for theta = theta_range
d_ki2 = tan(theta) * d_ki;
% 循环不同的 theta2 值
for theta2 = theta2_range
RR = rotationMatrix(S_AR, -theta2); %顺时针
d_tki = RR * d_ki2';
d_tki=d_tki';
r = d_tki + S_AR;
% 解出 t 的范围 判断 r 是否与集热器相交 地系坐标
a = r(1)^2 + r(2)^2;
b = 2 * (dot_AG(1) * r(1) + dot_AG(2) * r(2));
c = dot_AG(1)^2 + dot_AG(2)^2 - 49/4;
discriminant = b^2 - 4*a*c;
if discriminant >= 0
t1 = (-b + sqrt(discriminant)) / (2*a);
t2 = (-b - sqrt(discriminant)) / (2*a);
% 使用 t 的范围得到 z 的范围 地系坐标
z1 = dot_AG(3) + t1 * r(3);
z2 = dot_AG(3) + t2 * r(3);
z_min = min(z1, z2);
z_max = max(z1, z2);
if (z_min <= (h1+4) && z_max > (h1-4))%有交集
valid_points = valid_points + 1;
end
end
end
end
end
end
end
eta_sb = valid_point_count / (length(X) * length(Y));
Dot_Sum_light=length(theta_range)*length(theta2_range);%每个区域总光线数量
eta_trunc=valid_points/(valid_point_count*Dot_Sum_light);
eta=eta_cos*eta_sb*eta_at*eta_trunc*eta_ref;%计算光学效率 eta
Store(num,1)=eta;
Store(num,2)=eta_sb;
Store(num,3)=eta_cos;
Store(num,4)=eta_trunc;
end
%计算定日镜场的输出热功率 E_field
Ai=W_a*W_b;%定日镜镜面积
E_field=E_field+Ai*Store(num,1);%E_field=E_field+Ai*eta;
if num== 2836
disp(num);
end
end
E_field=E_field*DNI;
result=[result;mean(Store) E_field/(num*W_a*W_b)];
end
R=mean(result);
EEEE=[W_a W_b R(5)*W_a*W_b*num];
end

function [P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =shadow_range(i, h2) %i 是-S_i
i = i / norm(i);
% 圆柱参数
global h1;
R = 7/2; % 圆柱半径
O_top= [0 0 (h1 + R)]; % 圆柱上底面中心坐标
% 计算与入射光 i 垂直的直径方向
diameter_direction = cross(i, [0, 0, 1]);
diameter_direction = diameter_direction / norm(diameter_direction); % 单位化
% 上底面和下底面上与入射光垂直的直径的两个端点
P_top_1 = O_top + R * diameter_direction;
P_top_2 = O_top - R * diameter_direction;
P_bot_1 = P_top_1 - [0, 0, 2*R];
P_bot_2 = P_top_2 - [0, 0, 2*R];
% 计算在 z=h 时的交点
t_top_1 = (h2 - P_top_1(3)) / i(3);
P_top_1_intersect = P_top_1 + t_top_1 * i;
t_top_2 = (h2 - P_top_2(3)) / i(3);
P_top_2_intersect = P_top_2 + t_top_2 * i;
t_bot_1 = (h2 - P_bot_1(3)) / i(3);
P_bot_1_intersect = P_bot_1 + t_bot_1 * i;
t_bot_2 = (h2 - P_bot_2(3)) / i(3);
P_bot_2_intersect = P_bot_2 + t_bot_2 * i;
end
function isInside = pointInRectangle(x, y, P_top_1_intersect, P_top_2_intersect,P_bot_1_intersect, P_bot_2_intersect)
% 定义矩形的顶点
polygon = [P_top_1_intersect; P_top_2_intersect; P_bot_2_intersect; P_bot_1_intersect;P_top_1_intersect];
% 计算交点数
crossings = 0;
for i = 1:length(polygon)-1
    if ((polygon(i, 2) > y) ~= (polygon(i+1, 2) > y)) && ...
        (x < (polygon(i+1, 1) - polygon(i, 1)) * (y - polygon(i, 2)) / ...
        (polygon(i+1, 2) - polygon(i, 2)) + polygon(i, 1))
        crossings = crossings + 1;
    end
end
% 根据交点数判断点是否在矩形内
if mod(crossings, 2) == 1
    isInside = true;
else

    isInside = false;
end
end
function RR = rotationMatrix(axis, angle)
c = cos(angle);
s = sin(angle);
t = 1 - c;
x = axis(1);
y = axis(2);
z = axis(3);
RR = [
t*x*x + c, t*x*y - s*z, t*x*z + s*y;
t*x*y + s*z, t*y*y + c, t*y*z - s*x;
t*x*z - s*y, t*y*z + s*x, t*z*z + c
];
end


function data=data_func()% 这里求出的定日镜坐标是相对坐标
global W_a;
global L0;
R0 = 100;
delta_D=5+W_a+0.05;%给点容差%% !
coordinates = []; % 保存所有坐标
maxNumMirrors = 10000; % 估计的最大定日镜数
coordinates = NaN(maxNumMirrors, 2); % 预先分配空间
idx = 1; % 当前填充到的坐标索引
% 第一区域的五层
while true
R = R0 + delta_D;
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
% 重新计算 theta_D 的值
theta_D = 2*pi/ND;
% 计算 cos 和 sin 值
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
% 计算所有定日镜的坐标
x_values = R * cos_values;
y_values = R * sin_values;
% 更新坐标数组
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
i = 1;
while true
prev=sqrt((coordinates(idx-1,1) - coordinates(idx-2,1))^2 + (coordinates(idx-1,2)-coordinates(idx-2,2))^2);
R = R + sqrt(delta_D^2 - (prev/2)^2);
% 检查 R 的条件
if R > (350+L0)%%%%
break;
end
theta_D0 = atan2(coordinates(idx-ND,2), coordinates(idx-ND,1)) + theta_D/2;
cos_values = cos((0:ND-1) * theta_D + theta_D0); % 角平分线
sin_values = sin((0:ND-1) * theta_D + theta_D0); % 角平分线
x_values = R * cos_values;
y_values = R * sin_values;
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
distance_first_layer = sqrt((x_values(2)-x_values(1)^2+(y_values(2)-y_values(1))^2));
% 检查距离条件
distance = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) - y_values(1))^2);
if distance > 1.33 * distance_first_layer
break;
end
layer = layer + 1;
end
i = i + 1;
end
% 删除未使用的空间
coordinates = coordinates(~isnan(coordinates(:,1)), :);
%筛选定日镜坐标
% 检查约束条件
N = size(coordinates,1);
valid_indices = [];
for i = 1:N
x1 = coordinates(i,1);
y1 = coordinates(i,2);
if x1^2 + (y1-L0)^2 <= 350^2
valid_indices = [valid_indices; i];
end
end
data=[coordinates(valid_indices,1), coordinates(valid_indices,2)];
error2 = check_distance(data);
if ~error2
error('定日镜距离不符合');
end
data(:,3) = sqrt(data(:,1).^2 + data(:,2).^2);
data = sortrows(data, 3);
end
function result = check_distance(coordinates)% 检验点是否满足两个定日镜之间的距离要大于等于 delta_D
% 给定的参数
global W_a;
delta_D = 5 + W_a;
% 计算所有点之间的距离
distances = pdist(coordinates);
% 找到最小距离
min_distance = min(distances);
% 检查最小距离是否大于等于 delta_D
if min_distance >= delta_D
result = true;
else
    result = false;
end
end


